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Abstract: Bright point sources associated with extragalactic AGN and radio galaxies 
are an important foreground for low frequency radio experiments aimed at detecting the 
redshifted 21cm emission from neutral hydrogen during the epoch of reionization. The 
frequency dependence of the synthesized beam implies that the sidelobes of these sources 
will move across the field of view as a function of observing frequency, hence frustrating 
line-of-sight foreground subtraction techniques. We describe a method for subtracting 
these point sources from dirty maps produced by an instrument such as the MWA. This 
technique combines matched filters with an iterative centroiding scheme to locate and 
characterize point sources in the presence of a diffuse background. Simulations show that 
this technique can improve the dynamic range of EOR maps by 2-3 orders of magnitude. 
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1 Introduction 

Highly redshifted 21cm radiation emitted by neutral 
hydrogen gas during the Epoch of Reionization (EOR) 
contains a wealth of cosmological and astrophysical 
information. Several experiments aimed at detecting 
and characterizing this radiation are currently in progress 
or und er construction (MWAp] PAPEEJ^ LOFAE^ 
GMRT( [Pen et al.|2008| ). The EOR signal is also one 
of the main science targets of next generation radio 
facilities such as the SKj^ 

Observations of the Gunn-Peterson trough in the 
spectra of high-redshift quasars indicate that r eioniza- 
tion was essentially complete by redshift 2 = 6 ( Becker 
|et al.|[2001[ ). Meanwhile, the optical depth due to 
Thompson scattering observed for CMB photons by 
the WMAP satellite implies that for an instantaneous 
process Zreion = 10.9 ± 1.4 (Komatsu et al. 20091. 



Together, these observations support theoretical ex- 
pectation that reioniza tion was an extended process 
(Pritchard et al. 20091. The expected redshift range 
of reionization, z ~ 6 — 15, puts the 21cm line at 200-90 
MHz. Such low radio frequencies poses a number of ob- 
servational challenges. Terrestrial transmissions from 
radio, television and satellite communications are all 
prominent at or near this band. Additionally, refrac- 
tion of low frequency radio waves by the ionosphere in- 
troduces time- variable distortions into the observed ra- 
dio sky whi ch require continuous re-calibration ( Mitchef l 
et aH2008| ). 



-"^ http://www.mwatelescope.org/ 
^http://astro. berkeley.edu/ 'dbacker/eor/ 
^http://www. lofar.org/ 
^http://www. skatelescope.org/ 



Observations of the EOR signal are further compli- 
cated by astrophysical foregrounds. These foregrounds 
are dominated by galactic diffuse synchrotron emission 
(GDSE) and extragalactic AGN (point sources), with 
a smaller contribution from galactic free-free emission 



The GDSE has a steep spec- 
-2.5, [Rogers fc Bowmaiil|2008[ ) which 



(Shaver et al 

tral index (/3 

makes it hundreds of times brighter at 150MHz than in 

the 21cm rest-frame, while extragalactic poin t sources 

have a typical spectral index of j3 ~ —0.8 (Subrah- 



manyan 2002), with some evidence of flattening at 
lower frequencies (Cohen et al.|[2004 l. Interferometric 
observations are sensitive only to fluctuations in the 
sky brightness, and, taken together, these foregrounds 
are expected to have fluctuations at least three or- 
ders of magnitude greater than the EOR fluctuations 
( Santos et al. ,2005 ) . Initial observations of the dif- 



fuse foregrounds confirm t hese expectations ([Bernardi 
et al.|2009[|Ah et al. 12008) [Bernardi et al.|2010| ). 



Fortunately, the daunting task of isolating the com- 
paratively faint EOR signal is made tractable by the 
spectral smoothness of the foregrounds. The GDSE 
exhibits an essentially featureless power-law spectrum, 
while the combination of many point sources of vary- 
ing spectral indices can also be well-rep roduced as a 
smoothly- varying function of frequency (Wang et al. 
[2006). In contrast, observing the EOR signal at differ- 
ent frequencies corresponds to probing the relatively 
rapidly varying IGM density field across a substantial 
range of redshifts. Numerous authors have taken ad- 
vantage of this distinction to demonstrate the possibil- 
ity of subtracting foregrounds by fitting polynomials or 
other smoothly-varying functions al ong the observing 
frequency / line-of-sight dimension (Wang et al. [20061 
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Harker et al.||2009| |Bowman et a"L]|2009| [Gleser et aT] sky and is expect ed to be ~ 180K at 178 MHz ( jFurlari 



2008 Geil et al. 2008 \ 



etto et al. 



The effectiveness of these foreground removal strate- 
gies assumes the prior removal of the brightest point 
sources. Bright point sources cannot be effectively re- 
moved through line-of-sight subtraction because the 
necessarily incomplete itii-coverage of any real inter- 
ferometer inevitably creates a sidelobe pattern, evoca- 
tively termed "frizz" by Liu et al. (20091, across the 
plane of the sky. Changing the observing frequency 
changes the size of the synthesized beam and conse- 
quently moves this "frizz" across any given point on 
the sky. In this way, high angular-frequency struc- 
ture in the sidelobe pattern across the plane of the sky 
enters the observational frequency dimension, caus- 
ing "mode-mixing" in the line-of-sight power spectrum 
( Bowman et al.|2009 l. For this reason, most previous 
foreground removal studies assume there is some flux 
threshold Scut above which all point sources can be 
removed. Bowman et al. (2009) and Liu et al. (20091 



20061. For an isotropic power spectrurrri 



agree that Scut should be of order 10 — 100 mjy. How- 
ever, the method by which these bright sources are to 
be removed remains an open question. In this paper, 
we introduce a new technique for subtracting bright 
point sources from synthesis images produced by an 
instrument such as the MWA. 



this thermal noise is averaged over the thousands of 
baselines which sample each three-dimensional Fourier 
mode, leading to a formally highly significant detec- 
tion. 

The ability of the MWA to detect the EOR signal 
is in practice limited by the accuracy of the instrument 
calibration and foreground subtraction. As described 
in iMitchell et al. ( 2008 1 , the high data rate output by 



the MWA correlator precludes the long-term storage 
of the raw visibilities, necessitating real-time calibra- 
tion and imaging. The calibration cadence is set by the 
need to i) resolve temporal variation in the ionospheric 
refraction and ii) update the direction-dependent an- 
tenna response model as sources move across the sky. 
The MWA will complete a calibration and imaging cy- 
cle every 8 seconds. The calibrated images formed 
from these 8s integrations will then be coadded into 
stacks of ~ 8 minutes of observing time to further re- 
duce the data storage requirements. These co-added 8 
minute dirty maps form the basic data product from 
which the offline data analysis of the EOR experiment 
will be done. 



3 A Description of the Method 



2 The MWA EOR Experiment 3.1 Motivation 



The primary goal of the MWA EOR experiment is 
to measure the power spectrum of the cosmic neu- 
tral hydrogen density fleld. 21cm radiation from neu- 
tral hydrogen will be visible against the cosmic mi- 
crowave background (CMB) when the spin tempera- 
ture, Ts, deviates f r om the CM B temperature, Tcmb 
( Wouthuysenl|1952[ |Field||1959| ) with a pre dicted dif- 
ferential brightness temperature given by ( [Ciardi fc| 
Madau|2003"l ): 



(JTb = 26 mK XHi{l + 5) (l- 



TcM 

l + z 



0.02 



I 10 J Vn„J 



1/2 



(1) 



During reionization, Ts > Tcmb, and hence neutral 
hydrogen can be detected in emission. 

The detectability of the cosmological 21cm signal 
is fundamentally limited by the thermal noise of the 
radio telescope. For an interferometer, each visibility 
is subject to a thermal noise contribution given by (eg 
Thompson et al.|2001[ ) 



Vrn 



2kBTsi 



(2) 



where Tsys is the is the system temperature, Ac is 
the effective area of a single antenna, df is the chan- 
nel bandwidth, and r is the integration time. For the 
MWA, the system temperature is dominated by the 



Before describing our subtraction technique, we briefly 
review the treatment of point sources in a number 
of observational regimes which will motivate our ap- 
proach. 



3.1.1 Radio Astronomy 

The most well known and widely practiced method 
for removing point sources in radio astronomy is the 
CLEAN algorithm (Ho gbom|1974 | and its many vari- 
ants. CLEANing addresses one of the central difficul- 
ties in processing synthesis images; that it intrinsically 
involves deconvolution . The Cotton-Schawb CLEAN 
variant ( Schwab|1984 l, which subtracts the clean com- 
ponents from the ungridded visibilities and hence al- 
lows one to avoid aliasing and gridding errors, gen- 
erally produces the best results. Nonetheless, t he ob- 
tained dynamic r ange is usually limited to ~ 10^ ( Brigg j 
& Cornwell 1992 I. More sophisticated approac hes have 
been able to achieve dynamic ranges of ~ 10^ ( Cotton 



& Uson 12008 Voronkov & Wieringa 20041, however 



these techniques require real-time processing which would 
almost certainly exceed the MWA computing resources. 
Additionally, despite some theoretical efforts in quanti- 
fying the uncertainties associated with CLEANed im- 
ages (Schwarz 19781, practical experience has shown 
that CLEANing often alters image statistics and leaves 
spatially correlated residuals (ie 'stripes', [Cornwell et al.| 



^The neutral hydrogen density field is expected to un- 
dergo some cosmological evolution over the full 32 MHz 
MWA bandpass, but this does not affect the present dis- 
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1999 1 which would hkely corrupt measurements of the 



EOR signal. 

A second approach to subtracting radio point sources 



is 'Peeling' (Noordam 2004 van der Tol et al. 2007 



Interna et al.||2009l. Peeling is essentially the sequen- 



tial self-calibration and subtraction of the brightest 
sources in the field and is intended to overcome the 
calibration difficulties introduced by ionospheric re- 
fraction and direction-dependant gains which are in- 
evitable in wide-field low frequency radio observations. 
Peeling has the advantage of subtracting the point 
source model from the ungridded visibilities, similar to 
the Cotton-Schawb CLEAN variant. However, due to 
the need to update the calibration solution on timescales 
which are short compared to the timescale for iono- 
spheric variations, it is unlikely that more than ~ 
100 sources will be peeled by the MWA Real-Time 
Calibration, leaving hundreds of bright sources above 
the confusion limit. Additionally, the expected data 
rate from the MWA correlator (~ 19 Gb/s) precludes 
long-term storage of the raw visibilities, implying that 
post-processing to remove foregrounds will likely be 
restricted to the gridded visibilities / dirty maps. 

3.1.2 CMB 

Removal of point sources is also an important fore- 
ground subtraction step in processing of Cosmic Mi- 
crowave Background (CMB) temperature maps. Matched 
filters are used to increase the contrast between the 
point sourc es and the CMB or other di ffuse compo- 
nents (Togmark & de Oliveira-Costa 1998 1. CMB beams 



are sufficiently compact that it suffices to identify and 
mask out the bright pixels associated with the main 
lobe. 

3.1.3 Optical Astronomy 

Locating and measuring the flux of point sources are 
the most fundamental processes in optical astrometry 
and photometry. Images are usually convolved with 
the Point Spread Function (PSF) to maximize signal- 
to-noise. For reasonably oversampled images, it is rou- 
tine to obtain astrometric centroids on the order of 
10% of the pixel scale. ( |Pier et al.|2003 l. 



3.2 Subtraction Procedure 

Based upon the above considerations, we designed a 
procedure for subtracting bright point sources from 
synthesized radio images produced by an interferomet- 
ric instrument such as the MWA. As mentioned above, 
the very high MWA data rate and a finite storage ca- 
pacity imply that only time averaged dirty maps will 
be available for offline data analysis. Given a dirty 
map, our subtraction proceeds as follows: 

i) Filtering: We use a matched filter to opti- 
mize the signal-to-noise ratio and reduce contamina- 
tion in the flux and centroid estimates caused by dif- 
fuse sky emissio n (GDSE). As shown by Haehnelt & 
Tegmark (19961, the optimal linear filter, ip, is one 



non-pointlike components in the sky) is prominent. In 
Fourier space, this implies 



V^(k) oc f(k)/P(k) 



(3) 



We approximate r, the source profile, as a two- 
dimensional Gaussian fitted to the main lobe of the 
synthesized beam at the centre of the FOV. -P(k) is the 
power spectrum of the generalized noise. If there is no 
diffuse sky component, then P(k) = 1 and matched 
filtering is equivalent to convolving the dirty image 
with our Gaussian beam model. If there is a diffuse 
sky component, then we estimate P(k) iteratively from 
the data. Specifically, we first carry out our subtrac- 
tion procedure assuming that P(k) is constant. The 
result of this initial subtraction is a residual image with 
most of the point source flux removed. We then mea- 
sure P(k) from this residual image and fit the mea- 
sured values to a second order polynomial in log space 
(ie log(P(fc)) = a{log{k)f + b(log{k)) + c). We then 
repeat our subtraction procedure using this model of 
P(k) in Equation ^ In our simulations, we proceed 
to identify point sources as 5a local maxima in the 
filtered images. For future real data, it is likely that 
the locations of most bright sources would be known 
a priori from higher frequency observations. 

ii) Centroiding: We use the centroiding proce- 
dure developed for the Sloan Digital Sky Survey (SDSS, 
(York et al. 2000)) photometric pipeline to measure 



which upweights the source profile, while downweight- 
ing the scales at which the generalized noise (ie all 



the positions of detected sources in the filtered maps. 
This procedure uses Gaussian quartic interpolation to 
predict the centroid based on a local maximum and 
the eight surrounding pixels. There are two subtleties 
in this centroiding procedure: i) the main lobes are 
not actually Gaussian, or even symmetric. Hence, the 
measured centroids are biased and the bias is direc- 
tion dependent. We correct for this bias in the same 
manner as the SDSS; by simulating the beam at the 
inferred centroid position and determining the offset 
between the expected and measured positions ii) cen- 
troid and flux measurements of any given source are 
perturbed by the sidelobes of the other bright sources 
in the field. We address this problem through iter- 
ation; an initial estimate is made of the position and 
flux of each source neglecting the contribution of other 
sources, and this estimate is subsequently used to cor- 
rect for the sidelobe contribution at the position of ev- 
ery source. In principle, this process could be repeated 
until the parameters for each source converged, but in 
practice the generation of each source with sub-pixel 
positional accuracy is computationally intensive and 
the results presented in this work are based on a single 
such iteration. 

iii) Aperture Photometry: We estimate the flux 
of a source through the use of a circular aperture at 
the measured centroid position in the match filtered 
image. We correct for aperture bias by simultaneously 
measuring the aperture flux of a match flltered image 
of the beam reconstructed at the centroid position. 



4 Simulations 
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4.1 Description of Simulations 



We developed and tested our subtraction method through 
the use a simulation pipeline which combines MAPS, 
the MIT Array Preformance Simulator (Wayth et al, in 
preparation), with the MWA calibration and imaging 
system (MWA RTS, Mitchell 2008). MAPS is software 
package for simulating the observed visibilities gener- 
ated by a real interferometric array. In our case, it has 
been configured to simulate a 512 tile array with the 
antenna design and provisional layout of the MWA. 
Each simulation is a two second snapshot integration 
at an observational frequency of 178 MHz and over a 
channel bandwidth of 40 kHz. These simulations are 
computationally intensive as each integration contains 
over 10^ visibilites. Consequently, we make two im- 
portant simplifications at this stage. 

First, we ignore the effects of thermal noise. Each 
visibility should be subjected to a thermal noise contri- 
bution given by equation[2] However, this level of ther- 
mal noise is not appropriate for EOR foreground sub- 
traction as the foregrounds will be subtracted not from 
snapshots but rather from coadded maps which result 
from hundreds of hours of integration. Naively, we 
could introduce a longer integration time t to scale the 
thermal noise, but this does not account for the scaling 
of the sidelobe level which would result from the simul- 
taneous rotation synthesis. The correct scaling could 
be obtained by generating and coadding snapshots cor- 
responding to a ~ 6 hr exposure, however such a sim- 
ulation is beyond our present computational resources 
and is a subject for future work. 

Second, we assume that all MWA tiles have identi- 
cal response (gain). By extension, this assumption im- 
plies that the MWA beam in our simulations is direction- 
independant. On the other hand, our imaging proce- 
dure and point source subtraction technique do not 
make this assumption. Rather, they treat the data as 
in the general case of a direction-dependant beam. In 
future we intend to simulate a representative variation 
of tile gains, but in this work we only consider the sim- 
plest case of identical tiles. MAPS also has the ability 
to simulate the effect of the ionosphere on the observed 
visibilities for a given model sky, but in this work we 
neglect this effect. 

The simulated visibilities produced by MAPS are 
subsequently processed into synthesized images by the 
MWA RTS. The resulting SIN projection images have 
a pixel scale of 1.9 arcmin/pixel at the field centre. We 
cropped the central 512 by 512 pixels of each image 
(corresponding to a FOV of ~ 18 by 16 deg) and re- 
strict our analysis to this region. This same simulation 
procedure is also used to reconstruct the synthesized 
beam at any point in the field. The resulting synthe- 
sized beam at the field centre is illustrated in Figure 
[l] The sidelobe fluctuations seen throughout the held 
are ~ 1 % of the peak source fiux. 




Figure 1: The synthesized beam for a 512 element 
MWA. A logarithmic stretch has been appUed in 
order to emphasize the extended sidelobe struc- 
ture. 



5 Simple Tests 

5.1 Single Point Sources 

We begin by simulating 100 realizations of a randomly 
placed, isolated point source in an empty sky to mea- 
sure the accuracy of our centroiding algorithm in the 
ideal case. These point sources have a flux of 100 
Jy/beam, although the fiux normalization is not mean- 
ingful in this simple case. Applying our centroiding 
algorithm directly to the dirty maps yields an RMS 
centroiding error of 0.02 pixels. Measuring the cen- 
troid in the filtered images (P(k) = 1 in this case) 
yields an RMS centroiding error of 0.002 pixels. Fol- 
lowing correction for the centroiding bias described in 
Section |3.2[ we obtain an RMS centroiding error of 
10"'* pixels. These errors and all subsequent results 
are summarized in Table [T] 

We also measure the change in the dynamic range 
of our simulated images following point source sub- 
traction. The dynamic range is conventionally defined 
as the peak brightness divided by the RMS in empty 
parts of the image. We define empty regions as those 
pixels with counts less than twice the maximum value 
observed in our diffuse sky model (see below). In prac- 
tice, this threshold is about 3% of the maximum source 
flux. We are most interested in the relative dynamic 
range; the ratio of the dynamic range between the 
raw and subtracted images. Table IT] lists the dynamic 
range for the raw (unsubtracted) images, the Matched 
Filter subtracted images, and the Matched Filter and 
PSF bias corrected subtracted images. The listed val- 
ues are the logarithmic means of our 100 realizations. 
The relative dynamic range for the Match Filtered and 
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PSF bias corrected subtracted images is 2.0 x 10*. We 
illustrate the subtraction results for a single simulated 
source in Figure |2] Both panels in this figure show 
cross-sections through the brightest pixel in the field. 
The upper line in both panels is the raw image. We 
have plotted the absolute values to allow for a logarith- 
mic scale. The point source is clearly visible near the 
center of the field, as are the ~ 1% sidelobes. In the 
left panel, the lower line shows the residuals following 
subtraction from the Match Filtered image. The rel- 
ative dynamic range in this case is ~ 240. Note that 
the subtraction for this particular source happens to 
be substantially less accurate than the average of 100 
realizations reported in Table 11] (ie 7.9 x 10^). In the 
right panel, the lower line shows the residuals following 
subtraction from the Match Filtered image where the 
centroid has been corrected for PSF bias. The relative 
dynamic range in this case is ~ 7.6 x 10"^. In both 
cases, significant residuals remain at the position of 
the point source. The question of whether these pixels 
should be masked out or if they will be further reduced 
by continuum foreground subtraction is a subject for 
future work. 

5.2 Single Point Sources with a Dif- 
fuse Background 

Next, we add a diffuse sky component based upon de: 
Oliveira-Costa et al. (20081. We use a scaling appro- 
priate to an observational frequency of 178 MHz, and 
locate the center of our field at the planned center of 
the main MWA EOR field (RA 4h, Dec -26°). Hence, 
this diffuse component is identical for all of our simu- 
lations. Figure |3] illustrates the resulting diffuse emis- 
sion. As before, we simulate 100 realizations of a ran- 
domly placed, isolated point source which are added 
to this diffuse background. The Matched Filter in this 
case uses an empirical P(k) as described in Section 
|3.2| The RMS centroiding error following Matched 
Filtering and PSF bias correction is 5.5 x 10~* pixels. 
In order to be able to directly compare the results of 
these simulations to those performed with no diffuse 
component, we subtracted a simulated image of the 
diffuse sky without any point sources (ie the result of 
a perfect point source subtraction) from the residual 
image prior to calculating the dynamic range, so that 
the reported value is the dynamic range of the point 
source residuals only. To clarify, this subtraction of the 
diffuse sky is performed only after we have completed 
our point source subtraction procedure with the dif- 
fuse component present. The relative dynamic range 
for the Match Filtered and PSF bias corrected sub- 
tracted images is 2.5 x 10'^. 

5.3 Multiple Point Sources 

The centroid and flux estimates for each source are 
perturbed by the presence of all other sources in the 
held. We investigated this effect by repeating our pre- 
vious simulations with 10 randomly placed sources of 
equal (100 Jy/beam) flux in an empty sky. Equal flux 
sources correspond to the case of maximal mutual per- 
turbation. Each realization has a different random 




100 200 300 400 500 
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Figure 2: The results of point source subtraction 
for a single source. In the left panel, the upper 
line represents the raw image and the lower line 
represent the residuals following subtraction from 
the Match Filtered image. In the right panel, the 
upper line represents the raw image and the lower 
line represent the residuals following subtraction 
from the Match Filtered image with the PSF bias 
correction applied. See text for more details. 



distribution of sources within the FOV. The iterative 
correction for source centroids and fluxes described in 
Section |3.2| was applied. The RMS centroiding error 
following Matched Filtering and PSF bias correction is 
5.9 X 10~* pixels and the corresponding improvement 
in the dynamic range is 1.8 x 10^ . We note that the raw 
dynamic range in this case is predictably lower than in 
the single source case due to the increased sidelobe 
noise. 

5.4 Multiple Point Sources with a 
Diffuse Background 

We added our diffuse sky model to the realizations of 
10 randomly placed equal flux point sources. The RMS 
centroiding error following Matched Filtering and PSF 
bias correction is 5.5 x 10~* pixels and the correspond- 
ing improvement in the dynamic range is 1.9 x 10'^. 
Notably, the results are effectively the same with and 
without the diffuse component for these sufficiently 
bright sources. 

5.5 Comparison with Peeling 

We compared the results of our subtraction technique 
to the effectiveness of the MWA RTS peeling routine. 
As described above, peeling is performed on the un- 
gridded visibilites, theoretically allowing for a higher 
fldelity source subtraction. However, there are a num- 
ber practical considerations which limit the empiri- 
cal effectiveness of peeling. Peeling requires an in- 
put catalog of sources which are to be peeled. In 
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Figure 3: The diffuse sky component generated 
from de Oliveira-Costa et al. (2008). A linear 
stretch has been apphed in this case. 



ing algorithm will make use of previous solutions to 
improve the current peel. This improvement is not 
represented in our single snapshot simulations. 

For the ten sources in an empty sky simulations 
described in Section |5.3| the measured relative dy- 
namic range after peeling is 2.2 x fO'^. With the dif- 
fuse sky component added, the measured relative dy- 
namic range is 1.1 x 10''. We note that for these bright 
sources, the performance of our subtraction procedure 
from the dirty maps is comparable to peeling in the 
case of no diffuse background (1.8 x 10"^ compared 
2.2 X 10^) and actually performs better in the case of 
a diffuse background (1.9 x 10^ compared 1.1 x 10^). 
Although these examples do not constitute a definitive 
comparison of these techniques, they do illustrate that 
high-fidelity point source subtraction is not necessarily 
contingent on access to the ungridded visibilities. 



5.6 A More Realistic Source Distri- 
bution 

As a further illustration of our subtraction procedure, 
we created a simulations with and without diffuse back- 
ground together with 100 sources drawn from a popu- 
lation whose counts are given by 



a real survey, this catalog would be likely be gener- 
ated from high er frequency observations (eg Parkes 
Source Catalog, [Wright fc Otrupcek 19961. It would 
also be possible to be supplement the catalog with 
data from complementary low frequency facilities (ie 
GMRT) or, in a bootstrapping fashion, from earlier 
MWA observations. In any case, the input catalog 
would invariably contain flux and position errors. Ad- 
ditionally, the observed source positions are perturbed 
from their true positions by the ionosphere. To ac- 
count for these errors, the RTS peeling routine uses a 
local least-squares minimization to fine-tune the posi- 
tion and flux of peeled sources. For our simulations, 
the input catalogs used by the peeling procedure con- 
tain the true (simulated) positions and fluxes of the 
sources in the field, so that peeling with these values 
would produce a perfect subtraction. However, the 
local least squares minimization in this case actually 
serves to slightly perturb the point source parameters, 
leading to imperfect subtraction. Further, the ideal 
method of dealing with diffuse emission when peeling 
is uncertain. Peeling from the unweighted visibilities 
would produce flux errors of ~ 1% even for the bright- 
est sources. Presently, the MWA RTS implements a 
simple scheme in which baselines are weighted by a fac- 
tor of 1 — exp(— 6^/6o^), where b is the baseline length 
and 6o is a scaling parameter, both in wavelengths. We 
choose feo = 500, which effectively saturates the scaling 
and results in only the longest baselines making a sub- 
stantial contribution. Finally, due to the processing 
time constraints, MWA RTS peeling is non-iterative, 
meaning that the subtractions of the brightest sources 
are not updated to account for the subsequently in- 
ferred flux of the fainter sources. However, when pro- 
cessing a series of consecutive integrations, the peel- 



^(> Sjy) 



NoSjy-'^Jy 



(4) 



where we restricted the fluxes to be 1 — 100 Jy. For 
comparison, the 6C 151MHz counts ( Hales et al.|1988 l 
predict ~ 200 source (> IJy) across our fleld of view. 
The brightest simulated source happens to have a flux 
of 34.6 Jy. The sources were randomly placed on the 
sky, subject to a minimum separation of 10 pixels. We 
found that, unsurprisingly, our perturbative approach 
to centroiding produces signiflcant errors for close pairs 
of sources. For such sources a simultaneous source flt- 
ting procedure such as DAOPHOT ( Stetson 1987[ ) may 
need to be incorporated. Alternately, higher resolution 
radio data could be used to inform the centroid flts of 
blended sources. We have also ignored source cluster- 
ing, although it will be an important consideration in 
more realistic sky models ( |Di Matteo et al."[2002i ). Fig- 
ure |4] shows the resulting image for the case of sources 
with a diffuse background. 

For the case of 100 sources with no background, the 
RMS centroiding error following Matched Filtering is 
9.7 X 10~^ pixels and the corresponding improvement 
in the dynamic range is 190. The PSF Bias Correc- 
tion is effectively irrelevant due to the magnitude of 
the input centroid errors. Figure [5] shows the centroid- 
ing errors as a function of source flux. The centroiding 
errors are largely independent of source flux, imply- 
ing that the sidelobes of the brightest sources are ef- 
fectively being accounted for while measuring the cen- 
troids of the fainter sources. For the case of 100 sources 
with a diffuse background, the RMS centroiding er- 
ror following Matched Filtering is 0.04 pixels and the 
corresponding improvement in the dynamic range is 
80. Again, the PSF Bias Correction is effectively ir- 
relevant. The centroiding error is inversely correlated 
with source flux, as shown in Figure |6] The centroid 
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Figure 4: 100 Point sources together with a diffuse 
background as described in section [5. 6| 



fits of the fainter sources are significantly affected by 
the relatively bright diffuse background. 
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Figure 5: The centroiding errors (in pixels) as a 
function of flux for 100 sources with no diffuse 



background as described in Section 5.6 



5.6.1 Comparison to EOR Experiment Re- 
quirements 

Since we have access to the true fiux and position of 
each simulated source, we are able to examine the rela- 
tive dynamic range for each indivdual source. For each 
source, this quantity is simply the relative dynamic 
range of the residual image when the selected source 
has been subtracted with the fitted position and fiux. 
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Figure 6: The centroiding errors (in pixels) as a 
function of flux for 100 sources with a diffuse back- 
ground as described in Section [5^ 



as above, while all of the other sources have been per- 
fectly subtracted using the known true fluxes and posi- 
tions. In other words, while the above relative dynamic 
range values measure the factor by which the combined 
sidelobe level of all of the sources in the image has been 
reduced, the individual relative dynamic range mea- 
sures the factor by which the sidelobes of a particular 
source have been reduced. Figure [7| shows the individ- 
ual relative dynamic range as a function of source flux 
in the case of no diffuse background. It is evident that 
the dynamic range is correlated with source fiux; the 
brightest sources are least perturbed by the sidelobes 
of the other sources and their subtraction is relatively 
more accurate. This measurement also allows for a 
rough comparison to the values of Scut, the flux limit 
to which bright sources are required to be removed, 
as reported by Bowman et al. ( 2009 1 and Liu et al. 
(.2009). For example, suppose that a 1 Jy source was 
subtracted with a 1% flux-only error to produce a rel- 
ative dynamic range of 100. The residuals would then 
correspond to an unsubtracted 10 mjy source. Such a 
source is sufficiently faint that the foreground contin- 
uum subtraction procedure could tolerate its presence 
and still achieve the sensitivity required to detect the 
EOR signal. Although our subtraction procedure pro- 
duces both flux and centroid errors (as well as beam 
model errors in the case of calibration errors) , the rel- 
ative dynamic range provides an aggregate measure of 
the magnitude of the residuals. The two solid lines in 
Figure (tI indicate the relative dynamic range required 
to correspond to Scut levels of 10 and 100 mJy. Our 
subtraction exceeds the more optimistic 100 mJy level 
for all sources and approaches the 10 mJy level for 
most. Figure Is] shows the individual relative dynamic 
range as a function of source flux for the case with 
a diffuse background. As before, the dynamic range 
is correlated with source flux. In this case, nearly all 
sources exceed the 100 mJy cut level, but few exceed 
the 10 mJy level. 



Publications of the Astronomical Society of Australia 




10' 
Source Flux (Jy) 




10' 
Source Flux (Jy) 



Figure 7: The individual relative dynamic range 
as a function of flux for 100 sources with no dif- 
fuse background as described in Section [5^ The 
two solid lines represent the relative dynamic range 
corresponding to S^ut = 10 mJy (upper) and 100 
mjy (lower) 



Figure 8: The individual relative dynamic range 
as a function of flux for 100 sources with a dif- 
fuse background as described in Section [5^ The 
two solid lines represent the relative dynamic range 
corresponding to S^ut = 10 mJy (upper) and 100 
mJy (lower) 



5.6.2 Comparison with Peeling 

We also applied the RTS peeling, as in Section |5.5[ 
to our 100 source images. For the case of no diffuse 
background, the relative dynamic range after peeling 
is 4.6 X 10^. Peeling is more accurate for this exam- 
ple than in the 10 point source case, presumably be- 
cause the steep source counts imply that the brightest 
sources can be peeled first with relatively little side- 
lobe contamination. For the case of 100 sources with 
a diffuse background, the relative dynamic range after 
peeling is only 26. The residuals are dominated by a 
few intrinsically faint sources whose centroid estimates 
are stochastically degraded by the local least-squares 
minimization. This example should not be taken as 
demonstrating the limiting performance of RTS peel- 
ing; the peeling results could almost certainly be im- 
proved by some algorithmic fine-tuning. However, it 
does illustrate that a diffuse background will in all like- 
lihood significantly degrade the accuracy of peeling for 
fainter sources. 



6 Outstanding Issues 

Our results indicate that subtraction of point sources 
from dirty maps can substantially improve the sidelobe 
noise in wide field synthesis images. A simple estimate 
indicates that the accuracy of our subtraction attains 
or at least approaches the accuracy required for EOR 
detection as reported by previous foreground subtrac- 
tion studies. However, a number of open questions 
must be resolved before it will be possible to unam- 
biguously determine whether this approach can satisfy 
the requirements of an EOR experiment. In particular, 
outstanding issues include: 



Calibration Errors: The accuracy of the cal- 
ibration solution will be of vital importance to 
point source subtraction, both for determining 
source fluxes and for accurately reconstructing 
the sid elobe pattern of each source. |Datta et al.| 
( 2009 1 have recently investigated the calibration 



tolerances for an idealized peeling-like source sub- 
traction. A comprehensive study of the calibra- 
tion budget will need to account for a combi- 
nation of complex effects including time- varying 
calibrators, non-idealized dipole response, and 
the effects of the ionosphere. 

Multiple Frequencies: The MWA correlator 
will simultaneously generate visibilities for hun- 
dreds of frequency channels across at 32 MHz 
bandwidth. Recent work aimed at analyzing 



CMB data from the Planck satellite (Herranz 



et al. 2009 1 has demonstrated methods for ex- 



tending matched filters to point source detection 
in multi-frequency data. Application of these 
techniques to MWA data will effectively allow 
continuum fitting to be combined with the an- 
gular power spectrum to increase the contrast 
between point sources and the diffuse compo- 
nent. 

Extended Integrations: As mentioned in Sec- 
tion |4.1[ rotation synthesis will significantly alter 
(reduce) the sidelobes of all sources. Coadding 
images produced over an extended observing cam- 
paign will also probably require an additional 
level of regridding. It will also be necessary to 
model the time dependance of the synthesized 
beams. Presently, generating simulations which 
reproduce even a six hour integration is a signif- 
icant computational challenge. 
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Statistics of Residuals: We have used the dy- 
namic range as a simple metric for judging the 
effectiveness of our subtraction procedure. In re- 
ality, a more detailed understanding of the sub- 
traction residuals will be required to correctly 
ascertain the uncertainties in the measured EOR 
power spectrum. Bowman et al. ( 2009 1 have 



suggested the construction of statistical templates 
which could be used to distinguish foreground 
residuals from the EOR signal. 



Sky Model: The sky model of (de Oliveira- 



Costa et al. 2008) uses input data which limit 
its angular resolution to ~ 1 degree. Conse- 
quently, spatial power at small angular scales 
is underrepresented in our diffuse component. 
A method for introducing additional small scale 
power into this sky model is currently in devel- 
opment (Bowman, private communication). An 
alternative is to build up a generic sky model 
from known sources of emission as done by |Jelic| 
eFaLlpOOSl. 



• Out of Beam Sources: Bright sources which 
are located outside the primary beam can intro- 
duce significant sidelobe noise across the field of 
view. Our subtraction procedure, which relies 
upon the ability to locate the main lobe of the 
synthesized beam, may be inapplicable for such 
sources. Further, the direction-dependent an- 
tenna gains may be poorly constrained towards 
such sources. A satisfactory approach to deal- 
ing with out-of-beam sources remains an open 
question. 

7 Conclusions 

Bright point sources have previously been recognized 
as an important EOR foreground, but the method by 
which they should be removed has been unclear. In 
this work, we have presented a procedure for subtract- 
ing point source from radio interferometric synthesis 
images. We are able to increase the dynamic range 
of our simulated images by a factor of 2-3 orders of 
magnitude. The efficacy of this method relies in large 
part on the excellent uv coverage of the MWA. These 
results are comparable to the results of the RTS peel- 
ing, but are achieved from the dirty maps, alleviating 
the need for long-term storage of the ungridded visbil- 
ities. Of course, peeling is an essential element of the 
MWA calibration strategy, and hence these two tech- 
niques will be used in a complementary fashion; peel- 
ing will remove some number of the brightest sources 
in real-time, and a subtraction procedure such as we 
have described will then be used offline to subtract 
a larger population of fainter sources from the time- 
averaged dirty maps. Significantly larger image simu- 
lations, particularly over longer integrations and multi- 
ple observing frequencies, are required to further refine 
and validate this approach. Nonetheless, our initial es- 
timates indicate that this procedure will be able to re- 
move point sources with sufficient accuracy to satisfy 
the foreground subtraction requirements of the MWA 
EOR experiment. 
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Table 1: Point Source Subtraction Results 



Subtraction Step 



Ccntroid Error" Raw Dynamic Range Relative Dynamic Range 



Dirty Image 
Matched Filter 
PSF Bias Correction 



Single Point Source 
O02 
0.002 
10-'' 



No Background 
5.3 X 10^ 
4.2 X 10^ 
1.0 X 10^ 



1 
7.9 X 10^ 
2.0 X lO'' 



Single Point Source - Diffuse Background 



5.3 X 10^ 
1.6x10^ 
1.3 X 10^ 
No Background 



Dirty Image 
Matched Filter 
PSF Bias Correction 



0.02 
5.4x10-3 

5.5 X 10-4 

10 Point Sources 



Dirty Image 
Matched Filter 
RTS Peeling 

" In pixels. 1 pixel = 
~ 10- 



0.07 
0.04 



= 1.9 arcmin or 1 arcsec 
^ pixels. 



2.5 X 10^ 
2.0 xl03 



1 
2.9 X 10^ 
2.5 X 103 



Dirty Image 
Matched Filter 
PSF Bias Correction 
RTS Peeling 


0.02 1.7 X 10^ 
1.9x10-3 1.2x10'"^ 
5.9 X 10-4 3.1 X 10^ 




1 
6.9 X 102 
1.8 X 103 
2.2x103 


10 Point Sources - Diffuse Background 


Dirty Image 
Matched Filter 
PSF Bias Correction 
RTS Peeling 


0.02 1.8 X 10^ 
9.2x10-3 3.4x104 
5.5 X 10-4 3.3 X 10^ 




1 
1.9 X 102 
1.9 X 103 

1.1X103 




100 Point Sources - No Backgr 


ound 




Dirty Image 
Matched Filter 
RTS Peeling 


0.02 2.4 X 10^ 
9.7x10-3 4.4 xl03 




1 

1.9 X 10^ 
4.6x103 




100 Point Sources - Diffuse Back 


ground 





1 

80 
26 



